#include "WallImpingement.h"

WallImpingement::WallImpingement(Flag& pFlag, const shared_ptr<ModelSettings> modelSetting0,
	const shared_ptr<CaseSettings> caseSetting0, Particle& pIn,
	const shared_ptr<WallGeometry> wallGeometry0, const ImpingementParameter& ImpingePara_in)
{
	/*Section 1: Pre-setup*/
	double tol = 0.5 * pIn.diam;

	/*Section 2: Get the nearest distance, type*/
	/*Step 1: get the distance of particle to walls and outlet*/
	/*Step 2: compare the distance and type*/
	/*Step 3: Find the minimal one and the type*/
	bool outOfBC = false; // check if the particle is out of boundary;

	// bottom wall distance
	

	double shortBot = 1.0, shortTop = 1.0, shortOutlet = 1.0;

	//bottom
	double yBot = (*wallGeometry0).getBotY(pIn.pos[0]);
	if (pIn.pos[1] >= yBot)
	{
		shortBot = CalWallDist(pIn.pos, (*caseSetting0).coordinateVertex0, (*caseSetting0).coordinateVertex1);
	}
	else if (pIn.pos[1] < yBot)
	{
		outOfBC = true;
		shortBot = 0.0;
	}

	//top
	double yTop = (*wallGeometry0).getTopY(pIn.pos[0]);
	if (pIn.pos[1] <= yTop)
	{
		shortTop = CalWallDist(pIn.pos, (*caseSetting0).coordinateVertex3, (*caseSetting0).coordinateVertex2);
	}
	else if (pIn.pos[1] > yTop)
	{
		outOfBC = true;
		shortTop = 0.0;
	}
	
	// outlet
	if (pIn.pos[0] <= (*caseSetting0).coordinateVertex1[0])
	{
		shortOutlet = (*caseSetting0).coordinateVertex1[0] - pIn.pos[0];
	}
	else if (pIn.pos[0] > (*caseSetting0).coordinateVertex1[0])
	{
		outOfBC = true;
		shortOutlet = 0.0;
	}

	// Get the minimal distance
	if (!outOfBC)
	{
		vector<double> shortVec = { shortBot,shortTop,shortOutlet };
		ImpingePara.nearestBCFaceType = getMinPosition(shortVec);
		ImpingePara.shortDist = shortVec[ImpingePara.nearestBCFaceType];
	}
	else
	{
		ImpingePara.shortDist = 0.0;
		if (shortBot == 0.0)
		{
			ImpingePara.nearestBCFaceType = 0;
		}
		if (shortTop == 0.0)
		{
			ImpingePara.nearestBCFaceType = 1;
		}
		if (shortOutlet == 0.0)
		{
			ImpingePara.nearestBCFaceType = 2;
		}
	}

	/*Section 3: Get particle velocity reference to wall*/
	double directionP, angleP;
	vector<double> v0(3, 0.0);
	v0[0] = pIn.V[0]; v0[1] = pIn.V[1]; 
	vector<double> norm0(3, 0.0);
	if (ImpingePara.nearestBCFaceType == 0)
	{
		norm0 = (*wallGeometry0).normal_bot;
		vector<double> lineLeft= (*caseSetting0).coordinateVertex0;
		vector<double> lineRight= (*caseSetting0).coordinateVertex1;
		lineLeft.push_back(0.0);
		lineRight.push_back(0.0);

		ImpactAngleCorrect(norm0, lineLeft, lineRight, v0,
			directionP, angleP, pVelNormal2Wall, pVelTang2Wall);
	}
	else if (ImpingePara.nearestBCFaceType == 1)
	{
		norm0 = (*wallGeometry0).normal_top;
		vector<double> lineLeft = (*caseSetting0).coordinateVertex3;
		vector<double> lineRight = (*caseSetting0).coordinateVertex2;
		lineLeft.push_back(0.0);
		lineRight.push_back(0.0);

		ImpactAngleCorrect(norm0, lineLeft, lineRight, v0,
			directionP, angleP, pVelNormal2Wall, pVelTang2Wall);
	}
	else if (ImpingePara.nearestBCFaceType == 2)
	{
		pVelNormal2Wall = 0.0; pVelTang2Wall = 0.0;
	}

	/*Section 4: Check if particle is impingement and update particle flag and velocity is needed*/
	if (outOfBC) // out of domain 
	{
		if (ImpingePara.nearestBCFaceType == 2) //outlet and out of outlet
		{
			pFlag.reachOutlet = true;
			pFlag.UpdateFlag();
		}
		else
		{
			assert(false && "Particle is out of BC boundaries!");
		}
		
	}
	else
	{
		if (ImpingePara.nearestBCFaceType == 2) //closest to outlet
		{
			if (ImpingePara_in.shortDist > ImpingePara.shortDist && ImpingePara.shortDist <= tol)
			{
				pFlag.reachOutlet = true;
				pFlag.UpdateFlag();
			}
		}
		if (ImpingePara.nearestBCFaceType == 0 || ImpingePara.nearestBCFaceType == 1) // closest to walls
		{
            if (ImpingePara_in.shortDist > ImpingePara.shortDist && ImpingePara.shortDist <= tol)
			{
				pFlag.impingement = true;
				pFlag.UpdateFlag();
			}
			else
			{
				pFlag.impingement = false;
				pFlag.UpdateFlag();
			}
		}
			
	}

}

void WallImpingement::ImpactAngleCorrect(const vector<double>& norm0, const vector<double>& lineLeft, const vector<double>& lineRight, const vector<double>& vp,
	double& directionP, double& angleP, double& normalVel, double& tangVel)
{
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	// Method:
	// angle P is the angle of the line starting from left to right reference to origin, so atan(dy/dx);
	// direction is used to respect the wall normal direction. Wall normal direction may be not located in the same way 
	// that the original xy 2d coordinate systems used, it may be opposite as wall normal always point inside the domain
	// wall normal: positive means point out of wall surface. so positive normal velocity means particles are leaving 
	// wall boundary.
	// steps:
	// step 1: get the impaction angle. Positive means location at 1st and 3rd quadrant
	// step 2: get the direction: cross product of wall normall to n1 and n2ton1, cross product of velocity and n2ton1
	//         if they have same direction, the direction is positive; vice visa
	// step 3: use the angle, direction, velocity to calculate the normal and tangential velocity
	//         for normal, use absolute value of complex number results and then multiply direction
	//         for tangential, no correction needed.
	// reference: matlab source code global2wall_Vector, need correction
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

	// get the data points
	vector<double> n2_n1(3, 0.0); //n2 to n1
	n2_n1[0] = lineRight[0] - lineLeft[0];
	n2_n1[1] = lineRight[1] - lineLeft[1];
	n2_n1[2] = lineRight[2] - lineLeft[2];
	// get the impingement angle
	if (n2_n1[0] != 0)
	{
		angleP = atan(n2_n1[1] / n2_n1[0]); // no need transfer to degree
	}
	else if (n2_n1[0] == 0)
	{
		if (n2_n1[1] >= 0)
		{
			angleP = 0.5 * M_PI;
		}
		else
		{
			angleP = -0.5 * M_PI;
		}
	}

	// get normal to n1 and velocity to n1
	vector<double> normal_n1(3, 0.0), velocity_n1(3, 0.0);
	normal_n1[0] = norm0[0] - lineLeft[0];
	normal_n1[1] = norm0[1] - lineLeft[1];
	normal_n1[2] = norm0[2] - lineLeft[2];
	velocity_n1[0] = vp[0] - lineLeft[0];
	velocity_n1[1] = vp[1] - lineLeft[1];
	velocity_n1[2] = vp[2] - lineLeft[2];

	// get direction
	double cros_normall_n2 = normal_n1[0] * n2_n1[1] - normal_n1[1] * n2_n1[0];
	double cros_vel_n2 = velocity_n1[0] * n2_n1[1] - velocity_n1[1] * n2_n1[0];
	double cros_direct = cros_normall_n2 * cros_vel_n2;
	if (cros_direct > 0)
	{
		directionP = 1.0;
	}
	else if (cros_direct <= 0)
	{
		directionP = -1.0;
	}
	else
	{
		assert(false && "Unkonwn wall normal velocity processing issue, please double check!");
	}

	// get wall normal and tangential velocities
	WallVelocity(vp, directionP, angleP, normalVel, tangVel);
}
void WallImpingement::WallVelocity(const vector<double>& v0, const double& directionP, const double& angleP,
	double& normalVel, double& tangVel)
{
	double c0 = cos(angleP);
	double c1 = sin(-1.0 * angleP);
	tangVel = v0[0] * c0 - v0[1] * c1;
	// abs is used for the reason that wall normall may be not equal to the rotated y, but we have used direction
	// we know that if particle velocity location is same to wall normall, then the normal velocity should be positive
	// so we can used abosolute and then multiply direction.
	normalVel = abs((v0[0] * c1 + v0[1] * c0)) * directionP;
}